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We extend the gravitational self-force methodology to identify and compute new ©(/i) tidal invari¬ 
ants for a compact body of mass /r on a quasi-circular orbit about a black hole of mass M ^ p. In 
the octupolar sector we find seven new degrees of freedom, made up of 3-1-3 conservative/dissipative 
‘electric’ invariants and 3-1-1 ‘magnetic’ invariants, satisfying 1-1-1 and H-O trace conditions. After 
formulating for equatorial circular orbits on Kerr spacetime, we calculate explicitly for Schwarzschild 
spacetime. We employ both Lorenz gauge and Regge-Wheeler gauge numerical codes, and the func¬ 
tional series method of Mano, Suzuki and Takasugi. We present (i) highly-accurate numerical data 
and (ii) high-order analytical post-Newtonian expansions. We demonstrate consistency between 
numerical and analytic results, and prior work. We explore the application of these invariants in 
effective one-body models, and binary black hole initial-data formulations, and conclude with a 
discussion of future work. 


I. INTRODUCTION 

The prospect of ‘first light’ at gravitational wave detectors has spurred much work on the gravitational two-body 
problem in relativity. It is now a decade since the first (complete) simulations of binary black hole (BH) inspirals 
and mergers in numerical relativity (NR) [1]. Such simulations have revealed strong-held phenomenology, such as 
‘superkicks’ [2], and have provided template gravitational waveforms. Yet, it may be argued, numerical relativity has 
also highlighted the ‘unreasonable effectiveness’ of both post-Newtonian (PN) theory [3], and the Effective One-Body 
(EOB) model [4]. 

BH-BH binaries, and their waveforms, are described by parameters including the masses M, p, spins, orbital 
parameters {p, e), etc. The parameter space expands for BH-neutron star (NS) binaries - a key target for detection 
in 2016 [5] - as tidal interactions also play an important role [6, 7]. Semi-analytic models, such as the EOB model, 
allow for much hner-grained coverage of parameter space than would be possible with (computationally-expensive) 
NR simulations alone. In addition, effective models can bring physical insight [8-10]. For real-time data analysis it 
may be necessary to blend effective models with surrogate/emulator models [11, 12] and careful analysis of modelling 
uncertainties [13]. 

By design, the EOB model [14-18] incorporates under-determined functional relationships, which are ‘calibrated’ 
with PN expansions and numerical data. Recently, it was shown that invariant quantities computed via the Gravi¬ 
tational Self-Force (GSF) methodology [19-21] can be used for exactly this purpose [15, 22-26]. In fact, as the GSF 
methodology is designed to provide highly-accurate strong-field data in the extreme mass-ratio regime [27, 28], it 
provides complementary constraints to PN and NR approaches, which excel in the weak-field and comparable mass- 
ratio regimes, respectively [29]. Thus, new GSF data, nominally limited in scope to the extreme-mass ratio regime, 
p/M <C 1, may immediately be applied to enhance models of comparable-mass inspirals, required for data analysis 
at, e.g.. Advanced LIGO [5]. 

In recent years, a growing number of invariant quantities, associated with geodesic orbits in black hole space- 
times perturbed through linear order 0{p/M), have been extracted from GSF theory. For quasi-circular orbits on 
Schwarzschild, these include (i) the redshift invariant [30, 31], (ii) the shift in the innermost stable circular orbit [32], 
(iii) the periastron advance (of a mildly-eccentric orbit) [32, 33], (iv) the geodetic spin-precession invariant [25, 34, 35], 
(v) tidal eigenvalues [26, 36, 37], (vi) certain octupolar invariants [26, 37]. Recently, (i) has been computed for eccentric 
orbits [33, 38], and (i)-(ii) have been computed for equatorial quasi-circular orbits on Kerr spacetime [39]. 

In 2008, the GSF redshift invariant at 0{p/M) was compared against a post-Newtonian series at 3PN order 
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(i.e., 0{v^/c^)) [30]. Many further PN expansions have followed for invariants (i)-(vi) at very high PN orders 
[25, 26, 34-36, 40-42]. An ‘arms race’ between numerical (GSF) and analytical (PN) approaches has developed, 
enabling precise comparisons of high-order coefficients [35, 36, 40-42]. Such comparisons are invaluable in quality 
assurance, as they have been used to correct small errors in both GSF calculations [36] and PN expansions [35]. 
Furthermore, in the ‘experimental mathematics’ approach [40, 43], high-order PN coefficients may be extracted in 
closed (transcendental) form from exquisitely-precise numerical GSF calculations. 

The purpose of this paper is to classify and compute GSF invariants at ‘octupolar’ order, i.e., featuring three 
derivatives of the metric, or equivalently, first derivatives of the Riemann tensor. This sector has been previously 
considered by Johnson-McDaniel et al. [44] and Bini & Damour [26], among others [45-48]. Our intention is to 
provide a complementary analysis which extends recent GSF work on the dipolar (spin precession) and quadrupolar 
(tidal) sectors. We aim for completeness, by (i) seeking a complete basis of octupolar invariants, (ii) providing both 
numerical GSF data and high-order PN expansions at 0(/i/M). 

In outline, the route to obtaining invariants is straightforward: (1) in the GSF formulation, the motion of a small 
compact body is associated with a geodesic in a regularly-perturbed vacuum spacetime [49, 50]; (2) the electric tidal 
tensor Sab of the regularly-perturbed spacetime defines an orthonormal triad at each point on the geodesic; (3) the 
covariant derivative of the Riemann tensor Rabcd;e resolved in this triad gives a set of well-defined scalar quantities 
{Xi}l (4) the functional relationships yj(n), where is the circular-orbit frequency, are free of gauge ambiguities; (5) 
we define the ‘invariants’ Axi{Sl) to be the 0{p.) parts of the differences Xi(n,/i) — Xi(n,/r = 0). 

The article is organized as follows. In Sec. II, we introduce electric and magnetic tidal tensors of octupolar order; 
decompose in the ‘electric quadrupole’ triad; examine the ‘background’ {p = 0) quantities; and apply perturbation 
theory to derive invariant quantities through 0{p). In Sec. Ill we describe various computational approaches for 
obtaining the regular metric perturbation and its associated invariants. In Sec. IV we present our results, primarily 
in the form of tables of data and PN series. In Sec. V we outline two wider applications of our work. We conclude 
with a discussion of progress and future work in Sec. VI. 

Conventions: We set G = c = 1 and use the metric signature -1-2. In certain contexts where the meaning is clear we 
also adopt the convention that M = 1. General coordinate indices are denoted with Roman letters a, b^c,.. ., indices 
with respect to a triad are denoted with letters i,j, k,.. and the index 0 denotes projection onto the tangent vector. 
The coordinates (t, r, 0, 4>) denote general polar coordinates which, on the background Kerr spacetime, correspond to 
Boyer-Lindquist coordinates. Govariant derivatives are denoted using the semi-colon notation, e.g., ka-,b, with partial 
derivatives denoted with commas. Symmetrization and anti-symmetrization of indices is denoted with round and 
square brackets, () and [], respectively. 


II. FORMULATION 

A. Fundamentals 

1. Tidal tensors 


We begin by considering a circular-orbit geodesic in the equatorial plane of the regularly-perturbed vacuum Kerr 
spacetime gab with a tangent vector u“. From the Riemann tensor Rabcd (equal to the Weyl tensor Cabcd in vacuum) 
we can construct electric-type and magnetic-type ‘quadrupolar’ tensors. 


^ab — Racbd'^ j 

(2.1) 

Bab = RlcbdU^U^ 

(2.2) 

where Rai,ad = ^Sab^^'^f^d- We may also construct ‘octupolar’ tensors, 


^abc — Radbe\c^ ^ ; 

(2.3) 

Babe = Rldbe^uC 

(2.4) 

The quadrupolar tensors are symmetric {Sab = ^bo, Bab = Bba), transverse (SabU^ = 0 = BabU^), and traceless 
{B'^a = 0 in general, S°‘a = 0 in vacuum). Similarly, the octupolar tensors are symmetric, and traceless in the first 
two indices (as Rabcd = Rcdab and Rab = R'^acb ~ 0) in vacuum. By contracting the Bianchi identity (or its dual) 
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^abcd e + ^abde-c + ^abec d = 0 , we observe that the octupolar tensors are also traceless on the latter pair of indices, 
= 0 = in vacuum. Note however that the octupolar tensors are not symmetric in the latter pair of indices, 
in general. 


2. Tetrad components 

Let us now introduce an orthonormal tetrad {cg = u“,e“} on the worldline and define tetrad-resolved quantities in 
the obvious way, so that 

XiOj... — Xabc..Ti U Cj . . . , (2-5) 

where Xabc... is any tensor and i,j, k G {1, 2, 3}. The quadrupole components are spatial, Sqq = = 0 = = Bqq. 

The octupole components are spatial in first two indices, but not in general. We may then consider three types of 
octupolar terms, namely. 


EijQi and ( 2 - 6 ) 

and similarly for B. Here () and [] denotes the complete symmetrization and anti-symmetrization of indices. 

Note that £ij is real and symmetric, and thus its eigenvalues are real and its eigenvectors are orthogonal. Thus, we 
may select our triad to coincide with the electric-quadrupolar eigenbasis. In other words, we choose the triad in 
which £ij is diagonal. We choose 62 to be the vector orthogonal to the equatorial plane. 


3. Equatorial symmetry 

For circular equatorial orbits, the reflection-in-equatorial plane symmetry implies that many components are iden¬ 
tically zero. Namely, 

£12 = £23 = Bii — Hi 3 = B22 = B 33 = 0, (2.7) 

'£’112 = ^222 = '£^233 = ^123 = Oj (2-8) 

^111 = ^122 = ^133 = ^113 = ^223 = B 333 = 0, (2.9) 

with all permutations of these indices also zero. 


B. Classification of octupolar components 

Now we consider the three types of terms (2.6) separately, and show that Xijo and Xiljk] may be derived from 
dipolar and quadrupolar terms, whereas X(ijk) encode new information at octupolar order. 


I • £ijo and BijQ 

For circular orbits, we have = wcg, = 0 and = —a;e“ where w is the precession frequency with 

respect to proper time, defined by parallel transport observed from the electric eigenbasis (c.f. Ref. [34, 36]). As the 
quadrupolar eigenvalues are time-independent on circular orbits, the only non-trivial components are 

£izo = to (£ii — £ 33 ), Bi2o = —u;B23, ^ 230 =^^ 12 . ( 2 - 10 ) 

2- £i[jk] and 


By virtue of the the Bianchi identity, 

^a[bc] ~ 2 ^ Bdabc) ; 



( 2 . 11 ) 
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We now (i) project onto the tetrad, (ii) use that Bij = ^CjkiRoiki and £ij = recall that the 

tetrad components in the electric frame are constants for circular orbits. Thus all components are zero except 

^2[23] = '^l[31] = ^^^B23, 

^3[31] = 'S^2[12] = —-^i^Bi2, 

^1[12] = ^3[23] = 2^ {£ll — £ 33 ) I 

and permutations thereof. 


3- £(ijk) 

In general, £(ijk) and B(^ijk) each have ten components satisfying 3 trace conditions, i.e., seven independent compo¬ 
nents each. For circular orbits, 4 electric and 6 magnetic components are zero, respectively, leaving 6 and 4 non-trivial 
quantities satisfying 2 and 1 non-trivial gauge constraints. In other words, there are 10 quantities we may calculate 
(given below), satisfying 3 non-trivial trace conditions; thus, 7 new independent degrees of freedom at octupolar order. 


( 2 . 12 ) 

(2.13) 

(2.14) 


Other octupolar quantities may be written 
EOB theory [see Ref. [26], Eq. (DIO)] is K 3 + 

K3+ = + '^(^333) 


4- Additional invariants 

in terms of the set identified above. For example, a relevant quantity in 
= £(^abc)£^°'^‘^\ which may be expressed as 


3 (£2 


( 122 ) 


+ ^( 133 ) 


'^(^311) 


' '^(^322) 


^'^(130)’ 


(2.15) 


where £( 130 ) = \£i30- 


C. Circular orbits: Background quantities 


Below we give the values of the tidal quantities for circular equatorial geodesics on the unperturbed Kerr spacetime, 
i.e., for test-masses ()x = 0). Here, the orbital radius is rg and the orbital frequency is D = \fM A- a\fM') where 
a is the Kerr spin parameter and a > 0 (a < 0) for prograde (retrograde) orbits. 

The tangent vector and electric-eigenbasis triad have the components [51] 


[C/,0,0, HI/], (2.16a) 

[0,v^/ro,0,0], (2.16b) 

[0,0,l/ro,0], (2.16c) 

(2.16d) 

where U = y/M /(Hrg^^u), Ag = rg — 2Mrg -|- and 

= 1 - 3M/ro -f 2ay/Mlrl^‘^. (2.17) 

The spin precession rate is a; = y/Mro/vQ. 
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1. Quadrupolar components 


The (non-trivial) quadrupolar components are 


^11 

£22 

£33 

Bi 2 

We note that ^623 = 0 on the background. 


M 3 MAo 
r3 uVg 
2M 3 MAo 

'0 ^ '0 
M 

'0 

(1 - a/^/WFo) 

9/2 2 

Tg 


(2.18a) 

(2.18b) 

(2.18c) 

(2.18d) 


2. Octupolar components 


In the electric sector, 

- 9Mro - 12aA/Mn( + Iba^) (2.19) 

'^■( 122 ) = ^3rg — 2Mro — 16a\/Mro + 15a^ j (2.20) 

f(i33) = -A (srg - 7Mro + 4av^) , (2.21) 

where A = •\/AoM/(rgr;^). We note that £( 311 ) = £( 322 ) = '^^( 333 ) = 0 on the background. 

In the magnetic sector, 

6(211) = +C (4rg - SMrg + 7a^) -V (4r^ - TMrg + Sa^) , (2.22) 

6(222) = -C {3rl - 6Mro + 9a^) +V (Srg - 4Mro + 5a^) , (2.23) 

S( 233 ) = -C{ rl- 2Mro - 2a^) +V { rg - 3Mro) , (2.24) 


where C = 2M^/^/(rg^^^r;^) and V = 3aM/{rQv'^). Note that ^ 6 ( 123 ) = 0 on the background. 

The ‘derived’ quantities Sijo, £i[jk], etc., may be easily calculated using Eq. (2.10), Eq. (2.12)-(2.14) and Eq. (2.19)- 
(2.24). For example, in the Schwarzschild (a = 0) case, using Eq. (2.15) yields 

6M2(1 - 2M/ro)(15r2 - 46Mro + 42M2) 


K3+ = 


rg°(l — iM/r^Y 


(2.25) 


D. Circular orbits: Perturbation theory 

Here we seek expressions for the octupolar quantities in the regular perturbed spacetime gab + where gab is the 
Kerr metric in Boyer-Lindquist coordinates, and = 0(/i) is the ‘regular’ metric perturbation defined by Detweiler 
& Whiting [49]. We work to first order in the small mass fi, neglecting all terms at 0{tY), and noting that the regular 
perturbed spacetime is Ricci-flat. 

We take the standard two-step approach [30, 31, 36]. For a given geodesic quantity y (e.g. 5(iii)), we first compare 
X on a circular geodesic in the perturbed spacetime with y on a circular geodesic the background spacetime at the 
same coordinate radius r = Tq. Then, noting that tq itself varies under a gauge transformation at 0(fi), we apply a 
correction to compare y on two geodesics which share the same orbital frequency H. 

Following the convention of Ref. [36] , we use an ‘over-bar’ to denote ‘background’ quantities, so that barred quantities 
such as are assigned the same coordinate values as in Sec. IIC. We use 6 to denote the difference at 
i.e., 6 ef = e“ — ef. At 0{g), 6 may be applied as an operator with a Leibniz rule d{AB) = {6A)B + ASB. In general, 
such differences are gauge-dependent. To obtain an invariant difference, we introduce the ‘frequency-radius’ tq via 

Q = '/M + g'/M). 


(2.26) 
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Then, we write 


X(»’n)-X(rn) = Ax(ro)+0(/r^). (2.27) 

Here x(^o) has the same functional form as x on the background spacetime, with rg replaced by tq. As Ax is at 
0(/r), we may parameterize Ay using the ‘background’ radius rg, rather than rn, as rg — rn = O(^). Such 

relationships, Ax(rg), are invariant within the class of gauges in which the metric perturbation is helically-symmetric 
(implying that ^ = 0 at the relevant order). 


1. Perturbation of the tetrad 

We may write the variation of the tetrad legs in the following way. 


5u°‘ — + /^ggCg, 

3 

Sef = /3^ou‘" + ^ /3iie“. 
t=i 


(2.28a) 

(2.28b) 


with the coefficients (3ab = 0{fi) to be determined below. First, we note that /3gg and /Jgg may be found by recalling 
key relations previously established in GSF theory for equatorial circular orbits on Kerr spacetime [30, 52], namely, 


^rg + — 2 a\/ Mrg^ 


Su* 1 ^ fro I 2 , 2 

2 ^°°" 2 Vm ■ 

1 , 1/2 

—w = Tihoo — —— I ri — 2Mrg + a 
2 2M V ° 


)Fr, 

(2.29) 

\Fr. 

(2.30) 


Here /igg = and the radial component of the GSF is given by 


Fr = 

2 or 


(2.31) 


Hence we have /3gg = ^/igg and /Sgg = — 5 y Fr, where /igg = and Fr = fj- ^Fr = 0{fj-) is the (specific) ra¬ 
dial self-force. The diagonal coefficients fin follow from the normalization condition, (^gat + (e“ -I- Sef ) (cj -I- Se^) = 

6 ij. That is, fin = —\hn, where hn = (no summation implied). From orthogonality of legs 0 and 3, we 

have fiso = fios + hgg where /igg = By similar reasoning, /Igg = /igi and /3gi -I- /lig -I- hgg = 0. To eliminate 

the residual rotational freedom in the triad at we now impose the condition that the triad is aligned with the 

electric eigenbasis, i.e., that Sij is diagonal in the perturbed spacetime (so that, e.g., fig = 0). From this condition it 
follows that 


P 13 

hi 


{SR)io 30 — £llhi3 
fll - f33 

— (di?)l03g -I- fgg/lig 

fll — £33 


(2.32) 

(2.33) 


where 


(di?)i03g = -5i?a6cde>^e^nf 


( 2 . 34 ) 
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8. Perturbation of octupolar components 


Here we present results for the perturbation of the (symmetrized) octupolar components £{ijk) B(^ijk)- The 
electric components are 

(2.35a) 
(2.35b) 
(2.35c) 


<5f(iii) — (^Vi?)(ioioi) + ( hoo — 2^11 ) ^(111 ) + 2/3o372.ioi3i, 


<5f(122) — (<5Vi?)(10202) + (^^00 — 2^^^ ~ h22j £(122) + 2/3o3”^10232, 

'^'^■(133) = (^^^)(10303) + ^^00 ~ ~ ^33^ ^(133) + 2/3o3’^10333 + ^/^SO W (£ii — £33) , 

2 - - 

<^'^(113) = (10103) + g/^lOW (£11 — £33) + hl£{lll) + 2£i3£(i33), 

<^'£^(223) = (^V^)(20203) +^31^(122): 

<^■^^(333) = (<5Vi?)(30303) + 3/33i£(i33). 


where 


((5Vi?)(i0i0fc) = SRabcd-eU^u‘^e[‘'e'^je’'^ , 
%0j3k = Rabcd-eU^^efe^f‘e‘fe‘i}. 


The magnetic components are 


((5H)(2ii) — (^Vi?*)20101 + ( hoo — hii — 2^22 ) ^(211 ) + 2/?03’^20131j 


((5H)(222) — (^Vi?*)20202 + ( ^00 ~ 2^^^ ' ^(^22) + 2 £o3’^20232j 


where 


{6B)(233) — (^Vi?*)20303 + (^^00 ” 2^^^ ~ ^^33 j ^(233) + 2£o3^20333 + 
((5H)(i 23) = (^Vi?*)io203 +/5 i 3^(233) +/331^(211) + g/3l0^^12- 


(^Vi?*)(, 0 i 0 ) 0 ) = 


'R'i0j3k Rabcd;e^'^ 


-(b-d) -ia-c -e) 

'’3 6- 


(2.35d) 

(2.35e) 

(2.35f) 

(2.36) 

(2.37) 

(2.38a) 

(2.38b) 

(2.38c) 

(2.38d) 

(2.39) 

(2.40) 


3. Invariant relations 

As noted above, the coordinate radius of the orbit, r = ro, is not invariant under changes of gauge (i.e., coordinate 
changes at 0{fi)). On the other hand, the orbital frequency Q is invariant under helically-symmetric gauge transfor¬ 
mations. Following Eq. (2.27), we may therefore express the functional relationship between y S {£( 111 ),...} and H 
as follows, 

X(’’n) =x(’’n) + Ax(ro)+e>(/r^), (2.41) 

where ro is the frequency-radius defined in Eq. (2.26), and Ax = 0{fj,). Note that x(rn) denotes the ‘test-particle’ 
functions defined in Sec. IIC evaluated at tq. By definition, we have AH = 0. At 0{tL), 

or, making use of Eq. (2.29) and (2.30) and (5H/H = Su^/u^^ — 6u*/u*, 

= (2.43) 

In summary, Ax defined by Eq. (2.43), Eq. (2.35) and Eq. (2.38) are the invariant quantities which we will compute 
in the next sections. 


J^. Further quantities 


In Sec. IIB we wrote £ijo, £i{jk]i ^i[jk] in terms of quadrupolar tidal components, and the spin precession 

scalar w. If required, one may deduce the variation of these components by applying A as a Leibniz operator. For 
example, starting with Eq. (2.10), 


^'S’lso = Aw [Ell — E 33 ) + w (A£ii — A^aa). (2.44) 

Numerical data for the variation in the quadrupolar components Afu,..., AyB 2 i,... is given in Table I of Ref. [36]. 
We may compute Aw from the redshift and spin-precession invariants, AU and Atp, using 

Aw = ^AU - UnAi/j, (2.45) 

together with the data in Table III of Ref. [36] . 

Similarly, the variation Ait'a+, for example, can be found by applying A in this manner to Eq. (2.15). This can 
then be related to the quantity 5 k 3 +, whose post-Newtonian expansion was given to 7.5PN in Ref. [26]. Noting that 
K 3 + = r^JC 3 + and 


r = 


1 


Vl - 3M/ro 

we then have a relation between the first-order perturbations 

AK3+ 


1 + 2^00 


0[t?) 


K. 


= Sk 3 + + 2/1. 


00 - 


3-1- 


(2.46) 


(2.47) 


III. COMPUTATIONAL APPROACHES 


In this section we outline our methods for computing the octupolar invariants for a particle of mass /r on a circular 
orbit of radius tq in Schwarzschild geometry. Our approaches break into two broad catagories: (i) numerical integration 
of the linearized Einstein equation in either the Regge-Wheeler (RW) or Lorenz gauge and (ii) analytically solving 
the Regge-Wheeler field equations as a series of special functions via the Mano-Suzuki-Takasugi (MST) method. In 
both cases we decompose the linearized Einstein equation into tensor-harmonic and Fourier modes and solve for the 
resulting decoupled radial equation. In this section, and subsections that follow, I and m are the tensor-harmonic 
multipole indices, w is the mode frequency and we work with standard Schwarzschild coordinates {t^r,9,ip). For this 
section let us also define / = /(r) = I — 2M/r. We shall also use a subscript ‘0’ to denote a quantity evalulated at 
the particle. Finally, note that for a circular orbit about a Schwarzschild black hole the particle’s (specific) orbital 
energy and angular-momentum are given by 


rg — 2M 
3/ro[ro - 3M)’ 


Eo = ro 



(3.1) 


respectively. 

For calculations in the RW gauge there is a single ‘master’ radial function, 'i'lmui, to be solved for each tensor- 
harmonic and Fourier mode [53, 54]. For circular orbits the Fourier spectrum is discrete and given by w = Wm = mEl 
where El = is the azimuthal orbital frequency. Consequently, we label the RW master function with only Im 

subscripts hereafter. The full metric perturbation can be rebuilt from the and their derivatives [55]. For I > 2 

the ordinary differential equation that obeys takes the form 

= Sid{r - ro) + S 2 S'[r - rg), (3.2) 

where r* is the radial ‘tortoise’ coordinate given by dr^jdr = f~^ and U[r) is an effective potential. The effective 
potential used depends on whether the perturbation is odd or even parity. For the odd/even parity modes, equivalently 
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I + m = odd/even, the potential is given by 


Urir) = ^ [1(1 + 1) 


Utir) = 




/ 

j,2A2 


3M\ 18M^ 


2 A" ( A + 1 + — J + 


A + 


M 


(3.3) 

(3.4) 


respectively, where X = (I + 2)(l — l)/2 and A = A + SM/tq. The form of the source terms, Si, also differ for the even 
and odd sectors. Explicitly, the odd sector sources take the form [56] 


CO _ ‘^P'Tofo^O 


(3.5) 

(3.6) 


For the even sector, we have 


5? 

^ 2 ^ 


pq£o 

r-o/oA 


§|/2A - (A(A + l)r2 + 6AMro + 15M2) 
.^0 


{rlpq£o)Y*^{0,<j)), 


YLiOA) 


ApClfS (l-2)\ 

roSo (^ + 2 )! 


y:a 0 ,ci^), 


where we have defined the following expressions for convenience: 


(3.7) 

(3.8) 


X, j,{e,(j)) = sm0deYijn(6,(l)), 

Y, j>^(0, cj)) = (5^0 + sin0cos6>(9e + ^ ^^ sin^g) Yi^(6>, 0) 

_ ^ _ /o 

^ (A + 1)A- 

For the radiative modes I > 2,m ^ 0 we will construct homogeneous solutions to Eq. (3.2) either numerically or as 
a series of special functions, as outlined in the subsections below. For the static (I > 2,m = 0) modes, closed-form 
analytic solutions to the homogeneous RW equation are known. In the odd sector these can be written in terms of 
standard hypergeometric functions: 


(3.9) 

(3.10) 

(3.11) 


= x-^-\Fi(-l -2,-1 + 2,-21, x) (3.12) 

=x’' 2 Fi{l-l,l + 2,,2 + 2l,x). (3.13) 


where hereafter an overtilde denotes a homogeneous solution, a ‘-I-’ superscript denotes an outer solution (regular at 
spatial infinity, divergent at the horizon), a ’ denotes an inner solution (regular at the horizon, divergent at spatial 
infinity) and x = 2M/r. In practice, we need only solve the simpler odd sector field equations, and construct the even 
sector homogeneous solutions via the transformation [55] : 


= 


A -f A2 ± HuM 


, ,2 9M^(r-2M) 

^ + ^ + r2(rA + 3M) 


o± 




o± 

Im 


3Mf—^ 

dr 


(3.14) 


Note this equation holds for both static and radiative modes. 

We construct the inhomogeneous solutions to Eq. (3.2) via the standard Variation of Parameters method. As 
the source contains both a delta-function and the derivative of a delta-function the inhomogeneous solution and its 
radial derivative will both be discontinuous at the particle. Constructing the inhomogeneous solutions then becomes 
a ‘matching’ proceedure with the jump in the field and its derivative across the particle governed by coefficients Si 
and ^ 2 . Supressing even/odd notation, we define matching coefficients as follows: 


DL = 


Wlr, 


Si 

fo 


2 MS 2 


So 

\[/T- 


(3.15) 
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with the usual Wronskian defined as Wim = im^r^tm ~ Finally we construct the inhomogeneous 

solutions via 


^Ur) = Df^^Ur). (3.16) 

where the D^’s are constants for all values of r. 

To complete the metric perturbation in the RW gauge we use the I — 0 and I = 1 results of Zerilli [57]. Detweiler 
and Poisson expressed these contributions succinctly for circular orbits [58]. For the monopole and static dipole we 
have 


° ( - + 


1 


/ 


r tq — 2M 


0 (r-ro), 


ii=o X 


il—l.m—Q c\ r ' 2 , 

htv = - 2 M'Co sm ( 


r'^/rl r <ro 
1 /r r > Tq 


(3.17) 

(3.18) 

(3.19) 


where 0 is the Heaviside step function and all other components are zero. The Z = 1, m = 1 mode does not contribute 
to our gauge invariant quantities so we will not give the explicit expression for the non-zero hu , htr and h^-r components 
of this mode (but as a check we use the expressions, given as Eqs. (5.1)-(5.3) in Ref. [58], to check that the contribution 
from this mode to our invariants is identically zero). 

As well as working in the RW gauge we also make a computation in the Lorenz gauge. Our code is a Mathematica 
re-implentation of that presented by Akcay [59] and as such we refer the reader to that work for further details. 


A. Numerical computation of the retarded metric perturbation 

For our RW gauge calculuation, as discussed above, analytic solutions are known for the monopole, dipole and static 
(to = 0) modes. This only leaves the radiative modes (Z > 2 ,to 7 ^ 0) to be solved for numerically. Our numerical 
routines are implemented in Mathematica which allows us to go beyond machine precision in our calculation with ease. 
Given suitable boundary conditions near the black hole horizon and at a sufficiently large radius (we discuss below how 
we choose these radii in practice), we use Mathematica’s NDSolve routine to solve for the inner and outer solutions to 
the homogeneous Regge-Wheeler equation (3.2). Inhomogenous solutions are then constructed by imposing matching 
conditions of these functions at the location of the orbiting particle. 


1. Numerical boundary conditions 

In order to construct boundary conditions, we use an appropriate power law ansatz for rw at each of our 
boundaries, given by 


n+ 

,f,oo _ , itor. 

(3.20) 


(3.21) 


n—0 


Recursion relations for the series coefficients can be found by inserting our ansatz into the homogeneous RW equations, 
and choosing a maximum number of outer and inner terms rimax = n± gives us initial values for our fields at these 
boundaries. Inserting (3.20) and (3.21) into (3.2) for the odd sector, we find the following recursion relations: 


a„ = -*- [(Z(Z -I- 1) - n{n - l))a„_i + 2Muj{n - 3)(n - l)a„_ 2 ], 

Zn 

bn = , [{l{l + 1) + 2n(n - 1) - 3)a„_i - {n+ l)(n - 3)a„_2]. 

n[n — AiMuj) 


(3.22) 

(3.23) 
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As discussed above we do not need to solve the even sector field equations as we can transform from the simpler odd 
sector solutions using Eq. (3.14). 

For the inner homogeneous solutions, the convergence of the series (3.20) improves with increasing n_, and in 
practice we choose n_ = 35. The outer solutions require more care, as the boundary at infinity is an irregular singular 
point. Our expansion in Eq. (3.21) is an asymptotic series and, as such, the series is not strictly convergent in n for 
a fixed roo. The ansatz will initially show power law convergence with increasing n, but for sufficiently high n the 
series will begin to diverge. At this point it is no longer useful to add higher order terms. Note for a fixed max value 
n+ the series will still converge with increasing Too as expected. After analysing this behaviour, we take n_|_ = 100 to 
get the best boundary conditions. Given the boundary expansions as a function of rfj/ooi for fixed n±, we must then 
choose a location for our boundary sufficiently close to r* = ±oo to give the desired accuracy. Setting the final term 
in our ansatz to be of order 10“'^, where d is our desired number of significant figures, we choose as our boundaries: 

Too = (an+10‘')i/"+, (3.24) 

r//= 2M+(6„_10‘^)-fo”-. (3.25) 

The expansions (3.20) and (3.21) give the boundary conditions in terms of an arbitary overall amplitude, specified 
by Go and 6o- As we first construct homogeneous solutions we can set these amplitudes to any non-zero value, and in 
practice we choose qq = bo = 1. The amplitudes are then fixed by the matching procedure described above. 


S. Numerical algorithm 

In this section we briefly outline the steps we take in our numerical calculation in the Regge-Wheeler gauge. The 
Lorenz-gauge calculation follows a very similar set of steps [59]. 

• For each Im-mode with I > 2 solve the odd sector RW equation, even ii I + m = even. For the radiative 
modes {I > 2,m ^ 0) calculate boundary conditions for the homogeneous fields at ru/ao using Eqs. (3.20) 
and (3.21). Using the boundary conditions, numerically integrate the homogeneous field equation (3.2) from the 
boundaries to the particle’s orbit at r = rg. For the static modes {I > 2,m = 0) evaluate the static homogeneous 
solutions (3.12)-(3.13) at the particle. Store the values of the inner and outer homogeneous fields and their radial 
derivatives at rg. 

• For I > 2 and l+m = even transform from the odd sector homogeneous solutions to the even sector homogeneous 
solutions using Eq. (3.14). 

• For all modes with I > 2 construct the inhomogeneous solutions via Eq. (3.16). 

• For the I > 2 modes reconstruct the metric perturbation using the formula in, e.g.. Refs. [55]. 

• Complete the metric perturbation using the monopole and dipole solutions given in Eqs. (3.17)-(3.19). 

• Compute the retarded field l-mode (summed over m) contributions to the octupolar invariants using the formulae 
in Appendix A. 

• Construct the regularized Tmodes using the the standard mode-sum approach. The resulting contributions to 
the mode-sum accumulate rather slowly as l~^. 

• Numerically fit for the unknown higher-order regularization parameters and use these to increase the rate of 
convergence of the mode-sum with 1. This procedure is common in self-force calculations and is described in, 
e.g., Ref. [60]. 

• To get the final result sum over I and make the shift to the asymptotically flat gauge as discussed in Appendix B. 

For rg > 4M we set the maximum computed Amode to be Zmax = 80. This is sufficient to compute the octupolar 
invariants to high accuracy - see Sec. IV for details on the accuracy we obtain. For orbits with 3M < rg < 4M we 
find we need an increasing number of Z-modes to achieve good accuracy in the final results, and for orbits near the 
light-ring (located at rg = 3M) we set Zmax = 130 in our code - see Sec. IV C. 
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B. Post-Newtonian expansion 


The generation of analytic post-Newtonian expansions for the octupolar gauge invariants requires a calculation 
of the homogeneous solutions of the Regge-Wheeler equation for each £ mode. A general strategy for doing this 
was described in [61]. The calculation is broken into three sections: (i) the exact results of Zerilli give the £ = 0,1 
components of the metric - see Eqs. (3.17)-(3.19), (ii) certain ‘low-£’ values calculated using the series solutions of 
Mano, Suzuki and Takasugi and (iii) ‘high-£’ contributions using a post-Newtonian ansatz. In a recent paper [42] this 
approach was optimised and improved allowing extremely high PN orders to be computed, which otherwise are only 
accessible by experimental mathematics techniques [40, 41, 62, 63]. In the rest of this section we give a very brief 
overview of our technique and refer the reader to Ref. [42] for further details. 

The analytic MST homogeneous solutions are expressed using an infinite series of hypergeometric functions denoted 
which satisfies the required boundary conditions at the horizon, and a series of irregular confluent hypergeometric 
functions, X'^^, satisfying the boundary conditions as r» —>■ oo. Specifically we can write 

~ ntrans , _ 

^ ^trans -iwr, r- —nn 

so that, with oq = &o = 1 in Eqs. (3.20) and (3.21), we have the identihcation 

XZ = V (3-26) 

= (3-27) 

For the purposes of doing a PN expansion of the solutions for a particle on a circular orbit one finds two natural and 
related small parameters, the frequency w and the inverse of the radius, which are related by = ^J2GM/r^. A 
natural way to deal with this double expansion is to instead expand in 77 = 1 /c, and introduce two auxiliary variables 
Xi = GM/r, X^!"^ = wr, so that each instance of X^ and X^ must each come with an 77 ^ and are of the same order in 
the large-r limit. Expanding these solutions to a given PN order in this way amounts to truncating the X'"/“p infinite 
series at a finite order. However an in depth analysis of the series coefficients and the sometimes subtle behaviour of 
the hypergeometric functions reveals a structure that can be exploited to optimise this truncation order and fine tune 
the length of the expansion of each term in the series. 

A practical difficulty of this approach is that the MST series becomes increasingly large with higher 77 -order. For 
each PN order 7/ ~ W ^ ^2 gg.j. g^y PN, we need 20 77 powers. Significant further simplifications of these 

large series can be made rewriting the expansion as, for example, 


Y'in(MST) _ 






[l + ^2yl. +7744+77*^4 + ...] 


(3.28) 


where the Ai are strictly polynomials in Xi, A 2 . Since 2 X 1 X 2 ^/^ = 2GMu}, we see that t/j is r-independent allowing 
it to be essentially ignored as it will drop out with the wronskian during normalisation. We note that the purely even 
series in 77 includes some odd powers that appear at £-dependent powers, and with these we also get extra unaccounted 
for log terms. For instance, for £ = 2 the first odd term is at 77 ^^. 

As such, for a large-enough £ (dependent on the required expansion order), the homogeneous solutions become 
regular enough to instead use an ansatz of purely even powers as the solution of the RW equation. The details of this 
are described thoroughly in [42]. Once the homogeneous solutions of the Regge-Wheeler equation have been obtained, 
the even-parity solutions can be expressed using Eq. (3.14). This allows us to reconstruct the full metric perturbation, 
and from there our gauge invariant quantities, entirely from the Regge-Wheeler series solutions. 


IV. RESULTS 

In this section we present results for the octupolar invariants computed for circular orbits in a Schwarzschild 
background. More specifically, we present the six electric-type invariants defined in Eqs. (2.35a)-(2.35f), and the four 
magnetic-type invariants defined in Eqs. (2.38a)-(2.38d). In Sec. IV A we exhibit numerical data, and in Sec. IV B we 
supply post-Newtonian expansions. In Sec. IV C we examine the behaviour of the invariants in the approach to the 
light-ring at r = 3M. 
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rn/M 

Af (111) 



Af(122) 



Af(133) 



4 

-6.87640142 x lO”"^ 

5.634572704 x lO""^ 

1.24182872 x 1C 


5 

-1.3622429846 x 10“^ 

9.45418747546 x 10"® 

4.1682423703 x 10"® 

6 

-5.61141083923 x 10"® 

3.477232505498 

X 

10“® 

2.13417833373 x 10“® 

7 

-2.925643118454 

X 

10“^ 

1.701979164325 

X 

10“® 

1.223663954129 

X 

10"® 

8 

-1.710986615756 

X 

10“^ 

9.592475191788 

X 

10"^ 

7.517390965770 

X 

lO"'^ 

9 

-1.075500995896 

X 

10“^ 

5.890480037652 

X 

10"^ 

4.864529921306 

X 

10"^ 

10 

-7.120764484958 

X 

10-^ 

3.838876753995 

X 

10"^ 

3.281887730964 

X 

10"^ 

12 

-3.494517915911 

X 

10-^ 

1.847169590917 

X 

10"^ 

1.647348324994 

X 

10"^ 

14 

-1.913146810405 

X 

10-^ 

9.995537359206 

X 

10“® 

9.135930744843 

X 

10"® 

16 

-1.133949991793 

X 

10-^ 

5.879346730837 

X 

10“® 

5.460153187097 

X 

10"® 

18 

-7.141332604056 

X 

10“® 

3.682750321689 

X 

10“® 

3.458582282367 

X 

10"® 

20 

-4.718352028785 

X 

10“® 

2.423514347706 

X 

10“® 

2.294837681079 

X 

10"® 

30 

-9.514915883987 

X 

10“® 

4.835793521499 

X 

10“® 

4.679122362488 

X 

10"® 

40 

-3.040712519124 

X 

10“® 

1.538289606495 

X 

10“® 

1.502422912629 

X 

10"® 

50 

-1.252723439259 

X 

10“® 

6.321181929625 

X 

10 -" 

6.206052462967 

X 

10"'' 

60 

-6.064208487551 

X 

10 -" 

3.054930741569 

X 

10 -" 

3.009277745982 

X 

10 -" 

70 

-3.282027079848 

X 

10 -^ 

1.651475883984 

X 

10 -" 

1.630551195863 

X 

10 -" 

80 

-1.927657419028 

X 

10 -" 

9.691577786310 

X 

10"® 

9.584996403967 

X 

10"® 

90 

-1.205253640043 

X 

10 -" 

6.055682885677 

X 

10"® 

5.996853514753 

X 

10"® 

100 

-7.917190975864 

X 

10“® 

3.975890910078 

X 

10"® 

3.941300065787 

X 

10"® 

500 

-1.277421047615 x 10"^° 

6.392477321681 x 10"^® 

6.381733154472 x 10"“ 

1000 

-7.991970194046 x 10"^^ 

3.997657790884 x 10"^® 

3.994312403162 x 10"®^ 

5000 

-1.279743808249 x 10"^^ 

6.399252758926 x 10"^® 

6.398185323566 x 10"®® 


TABLE I. Sample numerical results for the conservative electric-type octupolar invariants. 


A. Numerical data 

We have employed two independent calculations in the Regge-Wheeler and Lorenz gauges: see Sec. Ill or Ref. [59] 
for details, respectively. Both codes are implemented in Mathematica, which allows us to go beyond machine precision. 
We find that the Regge-Wheeler and Lorenz gauge results for retarded field contribution to the invariants agree to 
around 22-24 significant figures. This high level of agreement, exemplified in Fig. 1, increases our confidence in the 
validity of the numerical calculation. 

In Table I we present sample numerical results for the three conservative electric-type invariants. Table II provides 
the results for the three dissipative electric-type invariants. As the computation of the latter does not involve a 
regularization step, the dissipative results are considerably more accurate than for the conservative results. Our 
numerical results for the three conservative and one dissipative magnetic-type invariants are presented in Table III. 


B. Post-Newtonian expansions 


As outlined in Sec. IIIB, we have made a post-Newtonian calculation of the octupolar invariants using a method 
which builds upon the work of Ref. [42]. This method allows us to take the expansions to very high order. Results 
at I5th post-Newtonian order are available in an online repository [64]. Here, for brevity, we truncate the displayed 
results at a relatively low order: 


Afdii) = + 8 y^ + 30y® - ( 

f 1604627 


4681^2 
IT 


6413231^2 

TT — 


1711 _ 

6 512 

159664 


'■)i + 


136099 

400 


6255 _2 
1024^ 


V 630 49152 " 105 "1 5 

AC — 4„,4 7„,5 o,,6 I /'1369 S)§r7 2\^7 . ( 

Ai.(i22) - 4l/ - 3?/ - 9?/ -I- (-g- 2048’^ >y + \ 

^ 132611239 _ 240298535 _2 173416 


18416 


log 2+ 4:^ logs- 


120960 


240298535 2 
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-7 + 
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53416 


35 


192 " ' 5 
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2048^ 

5 ' 

79832 
105 

7 


4096 

5 


log2 - 


1024 


logy)/-^-2/ 


265^2+ l(^^+2Cpi0g2 


525 
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5 logy)/ 

19/2^0(ylO) 
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86708 

315 


logy) / + 


109568 _„19/2 
525 


0{y 
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7200 
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rn/M A5(ii3) ^£( 223 ) A5(333) 

4 1.43018712098924 x 10“^ -6.81363125080514 x 10“^ -7.48823995908726 x 10“^ 

5 1.69051912392376 x 10“® -6.68228419170062 x 10“'‘ -1.02229070475370 x 10“® 

6 3.93615041880796 x lO^^ -1.40326772303052 x lO''^ -2.53288269577744 x lO^^ 

7 1.24851076918558 x lO^^ -4.16707182592131 x 10“® -8.31803586593453 x 10“® 

8 4.78575862364605 x 10“® -1.52593581419753 x 10“® -3.25982280944852 x 10“® 

9 2.09252624095044 x 10“® -6.45207930480653 x 10“® -1.44731831046978 x 10“® 

10 1.00921765694192 x 10“® -3.03317966765418 x 10“® -7.05899690176504 x 10“® 

12 2.90853534249746 x 10“® -8.42878520552755 x 10“^ -2.06565682194470 x 10“® 

14 1.02905687355232 x 10“® -2.90962875240999 x 10“^ -7.38093998311317 x 10“^ 

16 4.21307269886127 x 10“^ -1.17032511794287 x 10"^ -3.04274758091840 x 10“^ 

18 1.92451417988312 x 10“^ -5.27522804430828 x 10“® -1.39699137545229 x 10"^ 

20 9.57423553217574 x 10“® -2.59726822309717 x 10“® -6.97696730907857 x 10“® 

30 6.63075048503344 x 10“® -1.74627612621227 x 10“® -4.88447435882117 x 10“® 

40 1.00806706036123 x 10“® -2.61810134057605 x 10"“ -7.46256926303620 x 10"^° 

50 2.34720446527898 x 10"^“ -6.04700459200130 x 10"“ -1.74250400607885 x 10"^° 

60 7.14698583248068 x 10"“ -1.83158348544703 x 10"“ -5.31540234703365 x 10"^^ 

70 2.61736513652863 x 10"“ -6.68285229858544 x lO"^^ -1.94907990667008 x 10"^^ 

80 1.09692221205352 x 10"“ -2.79307929655166 x 10"^^ -8.17614282398351 x 10"^^ 

90 5.09523969822615 x 10"^^ -1.29465822904663 x 10"^^ -3.80058146917952 x 10"^^ 

100 2.56663032262918 x 10"^^ -6.51067248226772 x 10"^® -1.91556307440241 x 10"^^ 

500 7.32252735609857 x 10"^'^ -1.83582572460285 x 10~^'^ -5.48670163149571 x 10"^'^ 
1000 8.09167955607435 x 10"^® -2.02577720909791 x 10~^® -6.06590234697644 x 10"^® 
5000 2.31673725041822 x 10"®® -5.79347353907946 x 10"®'‘ -1.73738989651028 x 10"®® 


TABLE II. Sample numerical results for the dissipative electric-type octupolar invariants. 


rn /M _ AII(2ii) __ ^^{ 222 ) __ ^^( 233 ) __ ^^( 123 ) _ 

4 -6.148298254370 x 10"® 5.070286329453 x 10"® 1.078011924917 x 10“® 1.07801192491724 x 10"® 

5 -9.558323357929 x 10"® 7.670992694990 x 10"® 1.887330662938 x 10“® 1.88733066293848 x 10“® 

6 -3.155936380263 x 10"® 2.476758241817 x 10"® 6.791781384454 x 10“'‘ 6.79178138445361 x IQ-^ 

7 -1.397948966284 x 10"® 1.081923065789 x 10"® 3.160259004949 x 10“'‘ 3.16025900494923 x IQ-^ 

8 -7.234703923371 x 10"^ 5.550936330078 x lO""* 1.683767593293 x 10“^ 1.68376759329306 x IQ-^ 

9 -4.130973443372 x lO""' 3.151840931693 x lO^^ 9.791325116795 x 10“® 9.79132511679517 x 10“® 

10 -2.528015715619 x lO""' 1.921517184307 x lO^^ 6.064985313116 x 10“® 6.06498531311616 x 10“® 

12 -1.095551773983 x lO""' 8.288773695651 x 10"® 2.666744044177 x 10“® 2.66674404417714 x 10“® 

14 -5.444485917231 x 10"® 4.108569884562 x 10"® 1.335916032669 x 10“® 1.33591603266878 x 10“® 

16 -2.980264003325 x 10"® 2.245424872606 x 10"® 7.348391307187 x 10“® 7.34839130718667 x 10“® 

18 -1.753993694654 x 10"® 1.320134773357 x 10"® 4.338589212967 x 10“® 4.33858921296681 x 10“® 

20 -1.092394842331 x 10"® 8.215915676267 x 10"® 2.708032747040 x 10“® 2.70803274704046 x 10“® 

30 -1.770249199828 x 10"® 1.329249843091 x 10"® 4.409993567363 x 10'^ 4.40999356736253 x 10’^ 

40 -4.868817160862 x 10"'^ 3.653967138885 x lO"’’ 1.214850021976 x 10'^ 1.21485002197624 x 10’^ 

50 -1.788310960062 x 10"'^ 1.341778109443 x lO"’’ 4.465328506186 x 10“® 4.46532850618555 x 10“® 

60 -7.887212354333 x 10"® 5.917061308384 x 10"® 1.970151045949 x 10“® 1.97015104594935 x 10“® 

70 -3.946891664217 x 10"® 2.960771807198 x 10"® 9.861198570192 x 10“® 9.86119857019213 x 10“® 

80 -2.166444813367 x 10"® 1.625085708368 x 10"® 5.413591049991 x 10“® 5.41359104999132 x 10“® 

90 -1.276212037585 x 10"® 9.572758888543 x 10"® 3.189361487310 x 10“® 3.18936148731046 x 10“® 

100 -7.948907544564 x 10"® 5.962268324527 x 10"® 1.986639220037 x 10“® 1.98663922003664 x 10“® 

500 -5.716760192009 x 10“^® 4.287586670256 x 10"®® 1.429173521752 x 10~^® 1.42917352175245 x 10"^® 

1000 -2.528141980419 x 10“i® 1.896108307399 x 10"^® 6.320336730204 x 10~®^ 6.32033673020413 x 10"®^ 

5000 -1.809952182177 x 10“®® 1.357464188697 x 10“®® 4.524879934796 x 10“®'® 4.52487993479633 x 10“®® 


TABLE III. Sample numerical results for the magnetic-type octupolar invariants. 
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^'£’(113) 

^■^^(223) 

^■^^(333) 

^^(123) 

^^( 211 ) 


I / 62957089 394216079 _2 30557675496 ir.„o , 1458 i„„ q 152788 i .A ,,9 , 109 ^_,, 19/2 , 1 101 

+ l 17280 1179648 ^ 315 ' 35 ™ ^ + 7 315 ^°e,y)y + 525 

( 4 . 3 ) 

128 ^ 13/2 _ 128 ^ 15/2 ^ 5^^^8 _ mj^l 7/2 ^ 37 | 4 ^y 9 

+ 43^^(107554351 + 846720077 ^ - 258854407 - 51770880 log 2 - 12942720 log +0{y^°), ( 4 . 4 ) 

_ 32 13/2 _ 18 ,, 15/2 _ 128„„8 , 8276 „ 17/2 _ 15242 9 

5 y 5 y 5 « ^ 105 y 315 y 

- 43^(152535527 + 1693440077 ^ - 517708807 - 103541760 log 2 - 25885440 log + 0 (y^°), ( 4 . 5 ) 

_ 96 13/2 , 126 15/2 _ ^_,,8 , 38702 17/2 _ 3772 9 

5 y ^ 5 y 5 “y ^ 105 y 105 '^y 

- 45^(235966427 + 16934400772 - 517708807 - 103541760 log 2 - 25885440 log y)yi ®/2 o{y^°), ( 4 . 6 ) 

fi + f / + 

+ 116^(3475113 + 313600772 -9587207- 1917440 log 2-479360 logy)yi° - + 0{y^^), (4.7) 

-8y®/2 + Myll/2 _ 20yl3/2 + (-677 + f+ 772 )yl 6/2 

- 23^(246270016 - 20642025772 + 943718407 + 188743680 log 2 + 47185920 log y)yi ^/2 
- 1548^800 (417740314624 - 17848070625772 - 1319396966407 - 390644367360 log 2 

+ 1236197376001 og 3 - 65969848320 logy)yi ®/2 _ 2 ^ 6 ^^io (^(-^21/2^^ (-4 3^ 


AS(222) = 6y®/2 - 4y“/2 + ^yl3/2 + ^^2)^15/2 

+ 20^(234195584 - 19194125772 + 629145607 + 125829120 log 2 + 31457280 log y)yi^/2 
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-53747712000 log 3+ 51977256960 logy)yi®/2 54Z|4^^io (^(-^ 21 / 2 ^^ (4_4q-) 


where here y = M/tq. 

Figure 2 shows sample comparisons of our PN and numerical results. We observe that, as higher-order PN terms 
are included in the comparison, the agreement improves for all values of tq. For large orbital radii the comparison 
saturates at the level of our (smaller than machine precision) numerical round-off error. For strong-held orbits, the 
comparison allows us to estimate how well the PN series performs in this regime. At rg = lOM we typically hnd that 
the 15PN series recovers the hrst 7-8 signihcant digits of the numerical result. At the innermost stable circular orbit, 
at tq = 6 M, the 15PN series successfully recovers the Hrst 3-4 signihcant hgures. The excellent agreement we observe 
between our PN and numerical calculations gives us further conhdence in both sets of results. 


C. Behaviour near the light-ring 

With our numerical codes we can calculate the behaviour of the octupolar invariants as the orbit approaches the 
light-ring at rg = 3M. In general, the invariants will diverge as the light-ring is approached, and knowledge of the 
rate of divergence, along with our high-order PN results and our other numerical results, may be useful in performing 
global hts for the invariants across all orbital radii. Such hts hnd utility in FOB theory and already results for the 
redshift, spin precession and tidal invariants have been employed in FOB models [24-26]. In this section we discuss, 
and give results for, the rate of divergence of the invariants near the light-ring but stop short of making global hts for 
the invariants. 

The main challenge in computing conservative invariants near the light-ring is the late onset of convergence of 
the mode-sum in this regime (see Ref. [24] for a discussion of this behaviour). This necessitates computing a great 
deal more Zm-modes; typically we set Zmax = 130 for our calculations in this regime. By comparison, for orbits with 
rg = 4M we use Zmax = 80. Not only then do we need to numerically compute an additional 8085 Z?n-modes, on top of 
the 3239 modes required to reach /max = 80, but these higher /m-modes are more challenging to calculate numerically 
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FIG. 1. Comparison of numerical results computed in the RW and Lorenz gauges for a variety of conservative gauge invariant 
quantities, Axi, along a circular orbit at tq = lOM. We see 22-24 significant digits agreement in the individual tensor 1-modes 
of the retarded field. 
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FIG. 2. Comparison of our numerical and PN resnlts for (left) A£ii 22 and (right) AS 123 . For each invariant we plot the relative 
difference between the numerical data and successive truncations of the relevant PN series, i.e., in the legend ‘xPN’ means we 
are comparing against the PN series with all terms up to and including (relative) iPN order. As successive PN terms are added 
the agreement between the PN series and the numerical results improves. For the conservative invariants, such as A£i 22 , the 
agreement between the PN series and the numerical data saturates at a relative accuracy of 13-14 significant figures. For the 
dissipative invariants, such as Al?i 23 , the comparison saturates at 21-22 significant figures. This difference in accuracy in the 
numerical data stems from the requirement to regularize the conservative invariants whereas the dissipative invariants do not 
require regularization. 


owing to the stronger power-law growth near the particle for high I and the high mode frequency (and thus large 
number of oscillations that need to be resolved far from the particle) for high m-modes. These considerations mean 
that numerical calculations at radii near the light-ring are substantially more computationally expensive than our 
other numerical results. 

Our main results are presented in Fig. 3. We are able to infer the rate of divergence of five out of six of the electric- 
and magnetic-type invariants. Defining z = 1 — 3M/ro we find A£(iii) ^ —0.005892“®/^, A 5 (i 22 ) 0.00406z“®/^, 

^'^■( 133 ) ~ 0.0129z“®'^^, S( 2 ii) ~ —0.00392:“^ and ^ 8 ( 222 ) ~ 0.00392:“^ as z —)■ 0. For the remaining conservative 
invariant, A 8 ( 233 ), our current results are not sufficient to accurately determine the divergence rate, but we can say 
that the rate is subdominant to the other invariants. 


V. APPLICATIONS 


Here we briefly outline two possible applications of the results of Sec. IV : in informing FOB theory, and in refining 
initial data for binary black hole simulations with large mass ratios in the strong field. 
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FIG. 3. Divergence of the conservative octupolar invariants as the orbital radius approaches the light-ring. The electric-type 
invariants, Affm), A£(i 22 ), Af(i 33 ), diverge as where z = 1 — SM/tq. Two of the magnetic-type invariants, AI?( 2 ii) and 

AS( 222 ), are observed to diverge as z~^. We are unable to accurately deduce the rate of the divergence AS( 233 ) but we plot 
our numerical results to show that its rate of divergence is subdominant to the other invariants. 


A. Informing EOB theory 

In EOB theory, the dynamics of binary systems are reformulated in terms of the dynamics of a single “effective” 
body moving in a metric ds^ = —A{u; v)dt^ + B{u; i')df^+f'^ (dO^ -\- sin^ 6 d(fP) (non-spinning case), where A{u; v) and 
B(u; v) are smooth functions of inverse radius u = {M + ^)/f and symmetric mass ratio i/ = fj,M/(^ + M)'^. For tidal 
interactions, it was proposed in Ref. [65] that the metric function should take the form A = The 

latter terms are radial potentials associated with tidal deformations of bodies 1 and 2 , which may be decomposed into 
multipolar contributions, -|- A^^ ^ -|- A^ ^ from the electric quadrupole (A^^^), magnetic 

quadrupole (A^^”^), electric octupole (Ap'*’^), magnetic octupole sectors, respectively, etc. In Ref. [65] a 

relationship was established between the dynamically-sigmficant tidal functions A^^^^ and kinematically-mvariant 
functions J, (y) formed from the tidal tensors (see Eq. (6.11) in Ref. [26]). In the quadrupolar sector, the relevant 
invariants are 


_ c cab T _ >2 y 2 Q'b T _ c cbcco. 

q2 — ^ab^ 5 — f^ab^ t — ^ab^ 5 


(5.1) 


In the electric-octupolar sector, the relevant quantities are (see Appendix D of [26]) J 3 + = K 3 + + \J 2 +: where 
K 2 ,+ = and J 2 _,_ = In the magnetic-octupolar sector, analogous quantities may be formed. 

The 0{yL) part of these invariants may be easily deduced from our octupolar components A^m,.... For example, 
AK 3 +, obtained via Eq. (2.15), is related to 6 k 3 + by (2.47). 

Previously, Bini & Damour have given a PN expansion of ( 5 /^ 3 + to 7.5PN order (see Eq. (DIO) in Ref. [26]). With 
the results of Sec. IV, we are able to go a step further. First, in Table IV we give numerical data for Sk 3 + in the 
strong-held regime. The data indicates that 6 x 3 + has a local maximum somewhat within the innermost stable circular 
orbit. Second, in an online repository [64], we provide a higher-order PN expansion of Sk 3 +] below, we state the 
expansion at 8.5PN order (correcting a very minor error/typo in the y® term of (DIO) in Ref. [26]): 
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rn/M 

^3+ 

4 

-1.072402291940 

5 

-0.952268599881 

6 

-1.150905925689 

7 

-1.347913915585 

8 

-1.511472597166 

9 

-1.643850731891 

10 

-1.751437199028 

12 

-1.913557269058 

14 

-2.028682058336 

16 

-2.114109122984 

18 

-2.179795496907 

20 

-2.231771587180 

30 

-2.383995972376 
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-2.457665106706 
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-2.500976521370 
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-2.529455493583 

70 

-2.549596340465 
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-2.564588968234 
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-2.576181641423 
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-2.585412146067 

500 

-2.650685806947 
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-2.658693616512 

5000 

-2.665074853918 


TABLE IV. Sample numerical results for the Sk 3 + as defined in Eq. (2.47). 
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B. Informing initial data models 

How does a black hole move through and respond to an external environment? This question has been addressed 
by Manasse [66], and others [44, 45, 67-73], via the method of matched asymptotic expansions (MAE). In scenarios 
with two distinct length scales {M ^ /i), one may attempt to match ‘inner’ and ‘outer’ expansions across a suitable 
‘buffer’ zone {y r M) [74]. Indeed, this method was applied to derive the equations of motion underpinning 
the self-force approach [27]. Recently, much work has gone into improving initial data for simulations of binary black 
hole inspirals using MAEs [44, 46-48, 75, 76]. 

In a standard approach [67, 70, 71], the black hole is tidally distorted by ‘external multipole moments’: spatial, 
symmetric, tracefree (STF) tensors By, Eyfc, By^, etc., related to the Riemann tensor evaluated on the worldline 
in the regular perturbed spacetime. These STF tensors are essentially equivalent to our tetrad-resolved quantities; 
for example, Detweiler’s [70] STF moments are given by Ey = Sy, By = Hy, Ey^ = £(ijk) and Byj, = with 

the subtlety of the interchange of spatial indices 2 -o- 3. 

Johnson-McDaniel et al. [44] have applied the MAE method to ‘stitch’ two tidally-perturbed Schwarzschild black 
holes into an external PN metric. Implicit in Eqs. (Bla)-(Bld) of Ref. [44] is a PN expansion of (conservative) 
quadrupolar and octupolar tidal quantities. Restricting to 0{n), in our notation Eqs. (Bla)-(Bld) of [44] imply 

M3£-(ni)= 6 y'^ + 3y^ + -^{-8y^ + 8y^)+0{y^y^) (5.3) 

= -3y" - 42/' + ^ (4y" - 


( 5 . 4 ) 
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M3£:(i33) =-32/" + 2 /" + ^ (42/"-y2/') +0(2/®,/^^) (5.5) 

M^B^^ii) = 82 /®/^ + A (-8y9/2) + 0{y^^l\ /x") (5.6) 

M3S(222) = - 62 /®/^ + (62/®/2) + 0{y^^/\ y^) (5.7) 

^30(233) = -2y^/^ + ^ (22/'^/2) + 0(2/“/',/x^) (5.8) 


Note that here the 0{ijP) terms are leading-order terms in the Taylor expansion of the ‘background’ Schwarzschild 
results, and the 0{ij}) terms are consistent with the leading terms of our PN series in Sec. IVB. This reassuring 
consistency suggests that our 0{y/M) results may indeed be used to help improve initial data for large mass-ratio 
binaries in the latter stages of inspiral. 


VI. DISCUSSION AND CONCLUSION 

In the preceding sections we have pursued the line of enquiry of Refs. [26, 30, 34, 36], concerned with identifying 
and calculating 0{y) invariants for circular orbits, onwards into the octupolar sector. We identified 7 independent 
degrees of freedom in the octupolar sector, given by the (symmetrized) components of the derivative of the Riemann 
tensor as decomposed in the electric-quadrupole basis. A complete set of octupolar invariants for circular orbits is 
given by, e.g., A£(iii), Af(i 22 ), A,B( 2 ii), A,B( 222 ), A£i( 3 ii), Af( 322 ), A;B(i 23 ). Here, the first four are conservative and 
the latter three are dissipative in character. The remaining symmetrized components A£( 333 ), AB(^ 233 ) (conservative) 
and A£’( 333 ) (dissipative) follow from trace conditions. All additional octupolar components, A£yo, ABijo, ^£i[jk] ^^nd 
ABi^jk], may be written in terms of the previous-known quadrupolar tidal invariants ASu, A£ 22 , AS 12 , AS 23 [36], 
the spin-precession invariant A-tp [34] and the redshift invariant AU [30]. Accurate results for the latter quantities are 
provided in Tables I & III of Ref. [36] and PN series are given in Ref. [42]. In passing, we should note a relationship 
which was overlooked in Ref. [36]: AB 23 = —^12 Ax, where Ay is the dissipative invariant of Table I in Ref. [36]. 
Also, we should recall that the dissipative component of the self-force. Ft and are also invariants. Taken together, 
we believe we have now arrived at a complete characterization of all circular-orbit invariants in the regular perturbed 
spacetime through 0 (y), up to third-derivative order. 

Highly-accurate numerical results for all the octupolar invariants we identify are given in Tables I-IV. Our numerical 
calculation is performed using Mathematica and is made within the Regge-Wheeler gauge as described in Sec. III. In 
addition, as a cross-check on our results, we performed the same calculation in the Lorenz gauge using a Mathematica 
re-implementation of Ref. [59] - see Fig. 1 for an example of the excellent agreement we find between the two 
calculations. To complement our numerical results, we also calculate high-order post-Newtonian expansions for all 
the invariants. Our technique is briefly described in Sec. IIIB with the full details given in Ref. [42]. The lower-order 
PN expansions are given in Sec. IV B with the higher-order terms available online [64]. In Sec. V we explored two 
possible applications for the octupolar invariants. 

We can envisage several ways this work could be extended. First, the high-order post-Newtonian results and the 
strong-held numerical data could be combined to produce global semi-analytic hts for the various invariants. Here, 
knowledge of the behaviour at the light-ring (Sec. IV C) should prove useful. Similar hts for other invariants have 
already been applied to FOB models [24, 25, 37] and freshly-calibrated FOB models have been successfully compared 
against numerical relativity simulations [ 6 ]. Second, we note that in Sec. II we have, in fact, derived the form of 
the octupolar invariants for circular, equatorial orbits in a rotating black hole spacetime. Looking ahead, practical 
calculations on Kerr spacetime are needed. The redshift invariant has already been calculated for circular, equatorial 
orbits about a Kerr black hole [39, 52]. It seems a natural extension to extend other invariants, such as the ones 
we describe here, to the rotating scenario. We believe this should be pursued with both numerical and high-order 
post-Newtonian treatments. Third, a further natural extension is to consider invariants for non-circular orbits. This 
was recently explored by Akcay et al. [38] for the redshift invariant and we expect the calculation for other invariants 
to follow in time. Fourth, looking further into the future, invariants at second order in the mass ratio could be 
calculated. The necessary regularization procedure is now known [77-79] and the framework for making practical 
calculations is beginning to emerge [80, 81]. As with previous calculations, initial work will focus on the redshift 
invariant [82] but the calculation of other invariants will surely follow. 
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Appendix A: Gauge invariants in Schwarzschild coordinates 


In this Appendix we give explicit expressions for the perturbations to the octupolar invariants (as defined in 
Sec. IID 2) for the case of a circular orbit in Schwarzschild spacetime. Our expressions are written in terms of the 
components of hab and its partial derivatives in Schwarzschild coordinates, and are given by 
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Appendix B: Shift to asymptotically flat gauge 


In order to compare our results with PN theory it is necessary to work in an asymptotically flat gauge. In both the 
Lorenz and Zerilli gauges the tt-component of the metric perturbation does not vanish at spatial infinity and so we 
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make an 0{^) gauge transformation to correct for this [31]. For both gauges this correction can be made by adding 
= ^a-b + fb;a where = [—a{t + r* — r), 0, 0, 0] and a = fi/ \/ro(ro — 3M). Explicitly, this can be achieved by 
adding an extra term to the invariants, + 5^8 and + 5^8 where 
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